Commensurability and hysteretic evolution of vortex configurations in rotating optical 
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We present a theoretical study of vortices within a harmonically trapped Bose-Einstein condensate 
in a rotating optical lattice. Due to the competition between vortex-vortex interactions and pinning 
to the optical lattice we find a very complicated energy landscape, which leads to hysteretic evolution. 
The qualitative structure of the vortex configurations depends on the commensurability between 
the vortex density and the site density - with regular lattices when these are commensurate, and 
concentric rings when they are not. We model the imaging of these structures by calculating time- 
of-flight column densities. As in the absence of the optical lattice, the vortices are much more easily 
observed in a time-of-flight image than in-situ. 

PACS numbers: 37.10.Jk, 03.75.Lm 
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I. INTRODUCTION 

Atomic clouds in rotating optical lattices have gar- 
nered a large amount of interest from researchers in the 
fields of condensed matter physics, atomic physics, and 
quantum optics. An optical lattice simulates the periodic 
potential ubiquitous in solid state physics, while rotation 
probes the superfluid character of these cold atomic gases 
by driving the formation of quantized vortices. Here we 
explore the theory of vortices in an optical lattice. Specif- 
ically, we investigate the evolution of the vortex configu- 
rations that occur in the single-band tight-binding limit 
as the rotation rate is slowly varied. The energy land- 
scape has a complicated topography that leads to hys- 
teresis. The vortex configurations depend on commensu- 
rability of several different length scales. 

A uniform gas of atoms of mass m in an optical lat- 
tice rotating with frequency H is characterized by sev- 
eral important scales. Among these are the on-site in- 
teraction U, the lattice spacing d, the magnetic length 
i = \JhlmD,, and the particle density n, where h = hjl-K 
is Planck's constant. The behavior of the system changes 
when these various scales form different commensurate 
ratios. There are three well known examples of such 
commensurability effects, namely when cP /"kP is ratio- 
nal for a two dimensional noninteracting gas, when nd^ 
is an integer of a non-rotating gas in dimension D, or 
when 7rn£^ is rational for a two-dimensional lattice-free 
gas. The first example gives the Hofstadter butterfly 
single-particle spectrum the second the superfluid- 
Mott transition, and the third the fractional quantum 
Hall effect. Here we explore how the commensurability 
between € and d pla;Y?_9ut in an interacting superfluid, 
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We study the vortex configurations that emerge in a 
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harmonically trapped atomic cloud inside a rotating op- 
tical lattice in the single-band tight-binding limit. The 
resulting phenomenology is rich, as the vortex configu- 
rations depend on a number of factors, including: the 
vortex-vortex interaction, the vortex-pinning potential 
due to the optical lattice, the finite cloud size, and the 
past history of the cloud. Fast enough rotation of a 
uniform superfiuid results in the formation of an ar- 
ray of quantized vortex lines of cross-sectional density 
n^or, corresponding to a mean intervortex spacing of 

— Ij^pK. In an infinite system, these vortices ar- 
range themselves in a triangular lattice c onfig uration, but 
a finite cloud size produces distortions [3, [3j E^j [iBl- 
A co-rotating optical lattice introduces a vortex-pinning 
potential with minima at the optical lattice potential 
maxima (between the occupied sites). For commensu- 
rate vortex densities, this pinning can cause the lowest 
energy configuration to switch from a triangular vortex 
lattice, to one that shares the geometry of the optical 
lattice [13, [H, d, iH, m. In practice the vortices 
are insufficiently mobile to find the true ground state, 
and one typically sees some metastable configuration, for 
example with a number of domains separated by defects. 
We present a realistic simulation of these effects. 

We perform calculations in two-dimensions, modeling 
a harmonically trapped gas of bosons in a rotating square 
lattice. The two-dimensional case is convenient because 
we can then concentrate on the interaction between vor- 
tex cores in a single plane. This is also an experimentally 
relevant geometry, as the dimensionality of the system 
can be controlled by using an anisotropic harmonic po- 
tential, or optical lattice, where the hard trapping direc- 
tion is along the rotation axis of the optical lattice [23j . 
A recent experiment [l7j realized exactly this scenario 
by placing a rotating mask in the Fourier plane of a laser 
beam which forms an optical dipole trap. The mask con- 
tained three/four holes, producing a triangular/square 
lattice in the image plane, where the atoms were trapped. 
The lattice spacing was large due to the nature of their 
optics but can in principle be made small enough to ex- 
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plore the single-band tight-binding hmit that we inves- 
tigate. A similar setup, using a dual-axis acousto-optic 
deflector, promises to reach this limit in the near future 

We find hysteresis in our numerical algorithm, reflect- 
ing the complicated energy landscape for the vortex con- 
figurations, and discuss how similar hysteresis will appear 
in experiments. This landscape reflects the competition 
between vortex-vortex interactions and pinning to the 
optical lattice. In section II we describe our mean- field 
ansatz and numerical self-consistency routine. In sec- 
tion III we show how vortex configurations evolve from 
commensurate lattices to incommensurate ring-like struc- 
tures as the rotation rate is varied. In section IV we 
present the hystcrctic evolution of vortex configurations 
on spin-up and then spin-down. In section V we simulate 
the results of time-of-flight imaging of these systems, and 
in section VI we summarize our results. 

Previous work, focusing on the multi-band weak lat- 
tice limit, found vortex structures similar to those we see 
in our tight binding model [H, [l^ , but did not report 
on how these structures evolved as parameters were adi- 
abatically varied. Our discussion of the expansion of the 
rotating cloud initially in the single-band tight-binding 
regime is also novel. 



II. CALCULATION 

In the reference frame of the rotating optical lattice, 
our system is described by the rotating Bose-Hubbard 
hamiltonian 120, Sal : 



H 



E 



tiia]ai + h.c. ] + 



E 



u 



fii {fii - 1) - ^jn, 



where tij = 



t exp 



(1) 

is the hopping ma- 
The rotation vec- 



i /I' dr- A{r) 

trix element from site j to site i. 
tor potential, which gives rise to the Coriohs effect, is 
A{r) = [m/h) [vi X rjj = ttv (xy — yx), where v is the 

number of circulation quanta per optical-lattice site. The 
local chemical potential /ii = fio — m (lu'^ — 17^) in- 
cludes the centripetal potential. In these expressions, /io 



is the central chemical potential, uj is the trapping fre- 
quency, is the rotation speed, fi is the position of site 
i, m is the atomic mass, a\ (oi) is a bosonic creation 
(annihilation) operator, hi = a\di is the particle number 
operator for site i, and U is the particle-particle interac- 
tion strength. The connection between these parameters 
and the laser intensities are given by Jaksch et al. [2^ . 
Here, and in the rest of the paper, we use units where 
the lattice spacing is unity, and we operate exclusively at 
zero temperature. 

Both the superfluid and Mott insulator are well de- 
scribed by a spatially inhomogeneous Gutzwiller product 
ansatz 12511 . 
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where i is the site index, M is the total number of sites, 
\n)i is the n-particle occupation- number state at site i, 
and is the corresponding complex amplitude, with 
Tin l/nP = 1- This ansatz is more general than the more 
standard mean-field approximation described by the lat- 
tice Gross-Pitaevskii equation. In the limit where the 
latter works well, the two theories agree. The Gutzwiller 
ansatz has been used extensively to understand experi- 
mental results [13, nil , and is well suited for studying 
the vortex physics that we consider here. 

Using equation ([2|) as a variational ansatz^ we minimize 
the energy with respect to the {/^}. We then extract the 
density pi = X^™ '^l/nP ^^e condensate order param- 
eter = {d.i) = E„ at each site. The 
condensate density pi = |(ai)P is equal to the superfluid 
density in this model, and is generally not equal to the 
density. 

We use an iterative algorithm to determine the {/^} 
which (locally) minimize the energy. We use a square re- 
gion with L sites per side with hard boundary conditions. 
We find that we must take L much larger than the effec- 
tive trap radius so that our solutions do not depend on 
those boundaries. Typically we use 40 < L < 90. We cal- 
culate (Hrbh) using equation Minimizing (Hrbh) 
with respect to /" with the constraint — 1 = 

gives nonlinear eigenvalue equations, one for each site. 
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where the sum is over all nearest neighbors of site j, m 
is the particle-number index, Xj is a Lagrange multiplier. 



and Rjk = exp 



ijp' df-A{r} 



with i- 



1. We itera- 



tively solve these equations: first choosing a trial order- 



parameter field laj"'!, where aj = (dj), then updat- 

ing it by af = E„ V^f^U {{»?'"^}) Pn ([^T^]) ' 
where p is the iteration index. Similar calculations were 
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performed by Scarola and Das Sarma (29| to analyze the 
case where the single-particle Mott state is surrounded 
by a rotating superfluid ring. In the uniform case this 

as well as 
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algorithm has been used by Wu et. al. ^ 

Goldbaum and Mueller ^] and Oktel et. al. slISC 

Since equation ([3]) is highly nonlinear, we find that the 
solution that this iterative algorithm converges to is sen- 
sitive to the initial state we use. This feature allows us 
to see the hysteretic effects described in the text. Exper- 
iments will see similar hysteresis, but the precise details 
will differ from our calculations (for example the critical 
frequencies for vortex entry and egress will be somewhat 
modified). 

We systematically explore the phase space, varying the 
parameters in the hamiltonian. We simulate clouds with 
diameter from 30-60 sites, comparable to the sizes studied 
in experiments [1^ [131 . For the largest simulations we 
impose four-fold rotational symmetry, but introduce no 
constraints in the smaller simulations. 



III. COMMENSURABILITY AND PINNING 

We find a great variety of vortex patterns, including 
those resembling square vortex lattices. These are most 
stable at the rotation rates where they are commensu- 
rate with the underlying optical lattice. Commensurate 
Bravais lattices exist when 1/v is an integer, and com- 
mensurate square lattices when v — l/(v? + rn?), for 
integral n and m [13, [H, [ll] ■ Which vortex patterns ap- 
pear in a simulation, or in an experiment |3l| . depends 
on how the system is prepared. This hysteresis occurs be- 
cause the energy landscape has many deep gorges with 
near-degenerate energies, separated by high barriers. 

To illustrate this energy landscape, we fix t/C/ — 0.2 
and (Iq/U ~ 0.3, and study how the energy evolves as 
we vary the rotation speed. First, starting with the non- 
rotating ground state, we sequentially increase the rota- 
tion speed, using the previous wavefunction as a seed for 
our iterative algorithm. We adjust lo as we increase 
so that the cloud size, related to the Thomas-Fermi ra- 
dius, Rtf — \J m(3^-^ j remains effectively fixed. The 
energy as a function of rotation speed, shown in Fig- 
ure [TJa), has a series of sharp drops, corresponding to 
the entry of one or more vortices from outside of the 
cloud. At these rotation speeds the system jumps from 
one local minimum of the energy landscape to another. 

Figure [T] (b)-(g) shows the superfluid density and 
phases associated with the vortex patterns found dur- 
ing this adiabatic increase in rotation speed, where we 
impose four-fold rotation symmetry about the trap cen- 
ter. Subfigures (b) and (c) show a regular square vor- 
tex lattice seen near the commensurate ly — 1/(2 x 6^). 
Subfigures (d) and (e) show the vortex configuration at 
I' ^ 1/(2 X 3.76^) which is intermediate between the com- 
mensurate values v — 1/(2 x 3^) and v — 1/(2 x 4^). 
Rather than forming a square pattern, the vortex con- 
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FIG. 1: Adiabatic spin-up (color-online, one-column). 
Properties of cloud during adiabatic spin-up with parameters 
{tjU = 0.2, ^o/t/ = 0.3, Rtf = 15). (a) Energy vs rotation 
rate. Sharp drops indicate vortex formation. Energy scaled 
by on-site interaction parameter [/, and rotation rate quoted 
as V, the number of circulation quanta per optical-lattice site, 
(b), (d), (f) Superfluid density profile at parameters labeled 
in (a). Light-to-dark shading corresponds to low-to- high den- 
sity, and position is measured in lattice spacing. Light spots 
correspond to vortex cores. Red and green lines are guides 
to the eye. (c) , (e) , (g) Superfluid phase represented by Hue. 
Solid white circle denotes edge of cloud. Dashed lines are 
guides to the eye. Black circles denote vortex locations. In 
(b), (c) and (f), (g) rotation speed should favor a commen- 
surate square vortex lattice rotated by 7r/4 from the optical 
lattice axes, (d), (e) represents an incommensurate rotation 
speed. 
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figuration appears to be made of concentric rings. Such 
ring-like structures also occur for superfluids rotating in 
hard- walled cylindrical containers |32| , where boundaries 
play an important role. Despite this analogy, it appears 
that in the harmonic trap these circular configurations 
are not a consequence of the circular boundary. When 
we simulate an elliptical trap, we still find roughly cir- 
cular vortex configurations. As one increases v towards 

^ 1/(2 X 3^), a domain containing a square vortex lat- 
tice begins to grow in the center of the trap. As seen in 
subfigures (f) and (g), at commensurability the domain 
only occupies part of the cloud, even though one would 
expect that a uniform square lattice would be energeti- 
cally favorable. The inability of the system to find the 
expected lowest energy configuration during an adiabatic 
spin-up is indicative of the complicated energy landscape. 

The patterns that we find are largely determined by 
the symmetry of the instabilities by which vortices enter 
the system. For example, even when we do not impose a 
four-fold symmetry constraint this adiabatic spin-up ap- 
proach never produces square vortex lattices oriented at 
an angle other than it /A with respect to the optical lat- 
tice axes. On the other hand, we readily produce other 
commensurate vortex lattices by choosing the appropri- 
ate rotation speed and seeding our iterative algorithm 
with the corresponding phase pattern. We have veri- 
fied this approach with square vortex lattices oriented at 
various angles with respect to the optical lattice, taking 
^-1/2 ^ {4,5,6,7,8,9,10}, (5z/)-i/2 = {2,3,4,5,6} and 
(10i/)-V2 = {2,3,4}. 



IV. HYSTERESIS 

We further explore the history dependance of the vor- 
tex configurations by increasing, then decreasing Q. We 
do not impose a four- fold symmetry constraint, but take 
a smaller system with Rtf = 7. At any given f2, the 
energy shown in figure [2] (a) depends on the system's 
history. The increasing(blue)/decreasing(red) rotation 
curve has sharp energy drops signaling the introduc- 
tion/ejection of vortices to/from the system. The energy 
drops occur at different Q for spin-up and spin-down, in- 
dicating that the critical rotation speed for a vortex to 
enter or exit the system is different. Generically, there 
are more vortices in the system on spin-down than on 
spin-up. Depending on fi, one may find a lower energy 
state by increasing (subfigs. (b) and (c)) or by decreas- 
ing (subfigs. (d) and (e)) the rotation rate. As demon- 
strated by subfigs. (d) and (e), vortex configurations 
produced during spin-up typically have the four-fold ro- 
tational symmetry of the optical lattice, while the vor- 
tex configurations calculated during spin-down are more 
likely to break this symmetry. An experiment will dis- 
play the same qualitative features, but slightly different 
vortex configurations. 

When is changed more rapidly {i.e., the step-size is 
increased), we find more symmetry broken states than 
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FIG. 2: Hysteresis (color-online, one-column), (a) Energy 
versus rotation rate during increase (blue line) and decrease 
(red line) of v. Energy steps in the blue (red) curve correspond 
to nucleation (expulsion) of vortices, (b)-(e) Order parameter 
complex phase for parameters labeled in (a). Black circles 
are drawn around vortex cores, and white circles indicate the 
approximate extent of the gas. 

when we use small steps. The energy differences be- 
tween the symmetric and asymmetric states are ex- 
tremely small, so the energies shown in figure[2]are robust 
over a large range of sweep rates. 



V. TIME-OF-FLIGHT IMAGING 

In a previous paper Q we proposed detecting vortices 
in optical lattice systems through time-of-flight imag- 
ing [28, 33] , where at time t = one turns off the lattice 
and the harmonic trap, letting the cloud expand. Af- 
ter some fixed time t one then produces an absorption 
image of the cloud using a resonant laser beam. In a 
weakly-interacting gas, the density profile is related to 
the momentum distribution of the gas. Here we elabo- 
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rate on this argument, and show how the vortices will 
be observable in the time-of-flight (TOF) images. This 
complements other methods for extracting vortex prop- 
erties, such as Bragg spectroscopy "ss]. Recently Palmer, 
Klein, and Jaksch investigated time-of-flight expansion in 
the fractional quantum Hall limit [§]. 

We present a simple model where we neglect interac- 
tions during the time-of-flight. This approximation is 
quite good. First, the interactions between atoms from 
different sites can generically be neglected: by the time 
atoms from different sites overlap, the density is so low 
that they have negligable chance of scattering. Second, 
in the single-band tight-binding limit the kinetic energy 
from the zero-point motion of atoms on a single site 
should exceed the interaction energy, meaning that the 
trajectory of the atoms will only be slightly perturbed by 
the interactions. If one did include the effects of interac- 
tions during the expansion one would see a slight blur- 
ring of the interference patterns. This blurring comes 
from two effects: (1) atoms from sites with higher oc- 
cupation will be moving faster (the interaction energy is 
converted into kinetic energy), and (2) the interactions 
introduce phase shifts which depend on atom number. 

Within our approximation, calculating the density of 
the expanding cloud reduces to a series of single-particle 
problems. Taking the initial wavefunction to be given by 
equation after time t the wavefunction will be 
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where di{0) is the operator which annihilates the single- 
particle state in site i of the lattice. This operator evolves 
via the Heisenberg equation of motion, ihd^aj(t) = 

dj (t) , i?ticc where -fffrco is the Hamiltonian for non- 
interacting particles. This is equivalent to evolving the 
single-particle state annihilated by di(t) via the free 
Schrodinger equation. 

For this analysis we use the notation that f is a vector 
in the x — y plane, and z represents the coordinate in the 
perpendicular direction. We take the initial (Wannier) 
state at each site, (t>i (f, z), to be gaussian: 



(t>i (f, z) 



where A 



exp 



2A2 



and A 
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(5) 
and uj± 



being the small oscillation frequencies for each well. In 
the geometry we envision, uj±_ 3> tOosc- The wavefunc- 
tions at a time t after release of the trap are calculated 
by Fourier transforming (f, z) to momentum space, 
then time evolving under iJfrco and finally Fourier trans- 
forming back to position space to arrive at 0i {r,z,T) = 
(j>i (F, f{z,t), where the only thing we need to know 
about the transverse wavefunction f{z,t) is that it is nor- 
malized so / |/(z, t)pdz = I. The in- plane wavefunction 
is 



\7r (A2 -I- iht/m) 
and the column density of the expanding cloud is then 
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where Ui {rici) is the number of atoms (condensed 
atoms) initially at site i, and -(/'(?, z) is the bosonic field 
operator annihilating an atom at position (r, z). 

In the long-time limit where the expanded cloud is 
much larger than the initial cloud (i.e., Dj — ht/mX ^ 
Rtf)i this expression further simplifies, and one has 

n(r,t) = pir,t)[{N-N,) + \A{r,t)\^] (8) 

p(r,t) = {nDiy\-^''/^r (9) 

A(r,t) = ^a,e-^"-^^/^^\ (10) 

i 

where N and Nc are the total number of particles and 



condensed particles, respectively. The envelope, p{r,t), 
is a Gaussian, reflecting the Gaussian shape of the Wan- 
nier state. The incoherent contribution {N — Nc)p{r,t) 
has no additional structure. This is a consequence of the 
Gutzwiller approximation, which neglects short-range 
correlations. Adding these correlations would modify the 
shape of the background, but it will remain smooth. 

The interference term has the structure of the envelope 
p{r, t) multiplied by the modulus squared of the discrete 
Fourier transform of the superfluid order parameter. The 
discrete Fourier transform can be constructed by tak- 
ing the continuous Fourier transform of the product of 
a square array of delta- functions, and a smooth function 
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which interpolates the superfiuid order parameter. The 
resulting convolution produces of a series of Bragg peaks, 
each of which have an identical internal structure which is 
the Fourier transform of the interpolated superfiuid order 
parameter. The vortices will be visible in the structure 
of these peaks. 

Vortices in real-space lead to vortices in reciprocal 
space. This result is clearest for "lowest Landau level" 
vortex lattices [s^l which are expressible as an analytic 
function of z = a; + multiplied by a Gaussian of the 

I— 12 / 2 

form e~l^l , where w is a length scale which sets the 
cloud size. Aside from a scale factor and a rotation, the 
continuous Fourier transform of such a function is identi- 
cal to the original. More generally, the topological charge 
associated with the total number of vortices is conserved 
in the expansion process. 

Figure [3] displays simulated expansion images corre- 
sponding to the initial square vortex-lattice state shown 
in subfigures[T](a) and (b), where t/U = 0.2. In these im- 
ages, lighter colors correspond to smaller column densi- 
ties. Using rubidium-87 atoms on an optical lattice char- 
acterized by d = 410 nm and a hard axis lattice depth 
of 30 En, which are the experimental values in ref. [23|, 
we find that Vq = 5.7 En, which gives A = 84 nm. Sub- 
figures [3] (a)- (c) display absorption images: the left-hand 
side {x < 0) is the column density calculated with equa- 
tion ([7]), while the right-hand side {x > 0) shows this den- 
sity convolved with a 3 /xm wide Gaussian, representing 
the blurring from typical optics. Subfigure [3l^d) displays 
the long-time expansion limit column density calculated 
using equation (fTO|) . which only depends on the {/*i}'s 
and the ratio (A/d). 

Subfigure El^a) shows the in situ {t—0 ms) column 
density. At this stage the Wannier functions are tightly 
peaked about the lattice sites, resulting in a series of 
sharp density bumps. The heights of these bumps are 
slightly modulated due to the square vortex lattice: 
near the cores of the vortices there is a slight depletion 
of the density. Due to the small vortex size, none 
of this structure is seen once the image is convolved 
with the Gaussian. This demonstrates that a typical 
imaging setup would be unable to resolve the vortices. 
Subfigure [3ljb) shows the abosrption image after 2 ms 
time-of- flight. Several very broad Bragg peaks have 
developed, each showing a number of low density regions 
reflecting the square vortex lattice. Again, the vortex 



structure is lost upon convolution. Subfigure [Sjjc) shows 
an absorption image after 20 ms TOF where, even after 
convolution, the Fourier transform of the initial vortex 
pattern is clearly visible in each of the Bragg peaks. 
In their investigation of atoms in non-rotating optical 
lattices, Spielman et. al. [l^ allowed their atoms to 
expand for 20.1 ms before imaging. Subfigure [3t^d) is a 
column density calculated in the long-time limit using 
equation (fTO| . Aside from an overal scaling, and a slight 
rotation of the structure within each Bragg peak, the 
image after 20 ms is almost identical to the image seen 
in the long-time limit. 



VI. SUMMARY 

We have studied the vortex structures in a harmoni- 
cally trapped Bose gas in the presence of a rotating op- 
tical lattice. We discussed the hysteretic evolution of 
the vortex structures as the rotation rate is varied. This 
hysteresis is a signature of the complicated energy land- 
scape. We observed a tendency for the system to find 
regular lattices configurations when the vortex density is 
commensurate with the site density. At incommensurate 
vortex densities we instead see a circular vortex pattern 
which is robust against changing the aspect ratio of the 
trap. Finally we analyzed the time-of-flight expansion of 
one of these condensates. We find that the vortex pat- 
terns are readily observed in the structure of the resulting 
Bragg peaks. 
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FIG. 3: Simulation of column density seen in time-of flight absorption images. (Wide). Darker shading corresponds 
to higher column density. The initial state corresponds to the the 4x4 square vortex lattice pictured in subfigures [ljb),(c). 
(a)-(c) 0, 3, and 20 ms expansion. Left: column density profile calculated using equation (O; Right: column density convolved 
with a 3 /im wide Gaussian distribution to represent the finite resolution of a typical imaging system, (d) long-time scaling 
limit, where Dj = Kt/m\ S> Rtf- 
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